simon-5502-14-slides

What did we learn in biostats-2, module 02?

Three models

All use a continuous dependent (outcome) variables
All include multiple independent variables
Multiple linear regression (Week 02)
- All independent variables are continuous
Analysis of covariance (Week 03)
- Mix of continuous and categorical independent variables
Multi-factor analysis of variance (Week 04)
- All independent variables are categorical

Why three models?

Historical precedents
Different issues
- Multicollinearity
- Mediator variable
- Risk adjustment
- Moderator variable
- Interactions

If it seems arbitrary to have such different names, you are right. Some of this is historical. If you invent something, you get the right to name it, and different inventors for multiple linear regression, analysis of covariance, and multifactor analysis of variance just decided to use different names.

There are also issues that arise far more often in one model than the other. Multicollinearity, for example, is reserved primarily for multiple linear regression. You’ll hear discussion about multicollinearity this week as well as mediator variables. Risk adjustment is a topic for next week under analysis of covariance, and you will see discussion of moderator variables and interactions two weeks from now under multifactor analysis of variance. In theory, these issues are not exclusive to one model or another. It’s just easier to talk about them within a particular context.

The general linear model

Single model that unites all three models.
Use of indicator variables for categorical data
Not the same as general IZED linear model
- SAS: proc glm versus proc glim
- R: lm() versus glm()

It didn’t take long for researchers to discover a common linkage between multiple linear regression, analysis of covariance, and multifactor analysis of variance. In particular, you can treat categorical variables as if they were continuous through the use of indicator variables. You’ll see more discussion about this later.

One bad thing about statistics are all the terms that sound almost the same but mean something quite different. You may have noticed “moderator” and “mediator” from an earlier slide. I can never remember which is which. You’ll also learn about sensitivity and specificity in week 7 (Diagnostic testing) and these are also so easily confused.

But the worst, by far, are the general linear model and the generalized linear model. The IZED added to the end of “general” makes it an entirely different model.

You may see a bit about the generalized linear model in week 7 (logistic regression).

In any case, I have really preferred the use of a single model, the general linear model, in place of multiple linear regression, analysis of covariance, and multifactor analysis of variance. I am a “lumper” and not a “splitter”. For what it is worth, SPSS forces you to use the individual approaches at times instead of the single general linear model.

Arguments for the lm() function

formula = \(dependent-variable\) ~ \(independent-variables\)
- \(independent-variables\) can be numeric, factors, or strings
data =
subset =
na.action =
- na.fail
- na.omit
- na.exclude
other arguments

What did we learn in biostats-2, module 03?

What is a covariate?

Variable not of direct interest
- Relationship to outcome is already established
- Still must be accounted for
Examples
- Smoking in a cancer study
- Gestational in a neonatology study
A covariate can be continuous or categorical

What is covariate imbalance?

Difference in mean value between treatment and control group
- Often a problem in observational studies
- Sometimes a problem in randomized studies

Covariate imbalance is an issue in many research studies. It occurs in a comparison of an outcome between a treatment group and a control group (or maybe an exposure group and a control group). You want the treatment group to be identical to the control group in every way except for the treatment itself. But sometimes a covariate also differs between the treatment group and the control group. If that happens then the outcome could be influenced not by the treatment but by the covariate.

This is often a problem in observational studies.

It can happen in randomized studies at times as well. In theory, randomization will insure balance between the treatment and control group. Patients with large values of the covariate appear in the treatment and control group with equal likelihood. Patients with small values of the covariate appear in the treatment and control group with equal likelihood.

Sometimes, though, randomization doesn’t work. It relies on the law of large numbers and this doesn’t hold for some sample sizes. In particular, studies with less than 20 observations are fairly likely to see covariate imbalance, even with randomization. Even with larger sample sizes, sometimes you just get a bad batch of random numbers.

Why is covariate imbalance an issue?

Biased estimates
- Comparing apples to oranges
Harms study credibility

If a covariate is imbalanced, it can produce biased estimates. If younger patients are more likely to be in the treatment group, for example, and younger patients tend to have better outcomes, then you don’t know if the outcome variable is influenced by age instead of treatment. The variables are tangled up. This is comparable to the issue of collinearity in multiple linear regression.

Now the bias can go in either direction. Sometimes covariate imbalance will produce an artificial difference in the outcome between treatment and control. But it is also possible that a covariate imbalance masks a difference between the treatment and control.

You will find many times that the covariate imbalance does not produce any serious bias, but you still need to account for it. Failure to control for or to adjust for covariate imbalance will hurt the credibility of your study.

Examples of covariate imbalance

Age in a study of smoking and Down’s syndrome
Smoking in a study of artillery assignment and sperm count

A study of Down’s syndrome births had a covariate imbalance between women who smoked during pregnancy, the exposure group, and women who did not smoke during pregnancy, the control group. The average age in the exposure group was much lower than the average age in the control group.

I was involved in a study of lead exposure on male fertility. The exposed group were soldiers who worked on an artillery crew. These big guns could shoot missiles out at great speed, but when the missiles flew out, a cloud of lead dust washed back into the crew. Now, I should note the linkage between lead exposure and fertility is not well established. The control group were soldiers working in an office setting, far away from the big guns. It turns out that the proportion of smokers varied quite a bit between the exposed group and the control. When you are out in the field with explosives all around, smoking was totally banned. Now the artillery crew could still smoke while off duty, but the office workers had more opportunities to smoke on and off duty. This was in the 1990s before workplace bans on smoking were very common.

Now, I also have to admit that the link between smoking and male fertility is also not well established. But there was enough evidence to at least raise some concern about the covariate imbalance.

Covariate imbalance versus confounding

Covariate imbalance is simpler
Confounder definition relies on causation arguments

Preventing covariate imbalance

Randomization
Matching
Stratification

If you can, you should try to prevent covariate imbalance during the design of a research study. Randomization, if you can do it, is a great way to reduce the risk of covariate imbalance. It doesn’t always work, but it does work quite often. One of the nicest things about randomization is that it prevents covariate imbalance among those variables that you have measured, but it also prevents covariate imbalance among covariates that you didn’t measure, either because it was difficult or impossible to measure them.

Matching is another research strategy that helps to prevent covariate imbalance. For every treatment subject of a certain age, find a control subject with a closely matched age. Make sure that males are matched with other males and females are matched with other females.

You might choose other variables to match on. Just make sure that the covariates that you match on are important influences on the outcome. And don’t choose so many variables to match on that you have difficulty finding good matches.

Stratification also can help. Divide your patients into broad categories such as young, middle-aged, and older. Then randomly assign to treatment and control within these broad categories.

Adjusting for covariate imbalance

Propensity score models
Analysis of covariance

Variables cannot be in the causal pathway

Fixed at time of randomization
Temporally preceding exposure
Example: bottles given during a breast feeding study

Covariates must be fixed and in place at the time when you flip the coin to choose whether they get into the treatment group or the control group.

If you are studying an exposure, then the covariate must be a variable that precedes the exposure in time.

Variables that are intermediate between the treatment/exposure and the assessment of the outcome are said to be in the causal pathway.

If you adjust for variables in the causal pathway, that may reduce or even eliminate the effect of your treatment or exposure.

I saw an example of this in a study I helped with. It was an examination of breast feeding patterns in pre-term infants. It is difficult to maintain breast feeding in a pre-term infant because the mother goes home from the hospital before the baby does. A rule of thumb is that the number of weeks that a pre-term baby has to stay in the hospital is roughly equal to the number of weeks early that the baby arrived.

In this study, mothers of newborn infants were randomly assigned to a treatment or control group. In both groups, mothers were encouraged to breast feed when they were in the hospital visiting their baby. They were given breat pumps to collect milk when they were at home. The difference was that in the control group, infant were fed by bottle when the mothers were not around. In the treatment group, the infants were fed through an nasogastral (ng) tube. The intervention was designed to avoid having the infants becoming habituated to the artificial nipple of the bottle and then having trouble latching onto the mother’s nipple.

The intervention was quite successful. There were a few minor mistakes made during the experiement, though. Sometimes an infant in the treatment group was given a few bottles at the hospital instead of being fed exclusively through the ng tube. That’s not surprising and not too much of a cause for concern.

The researchers did measure the number of bottles given in the birth hospital in both the treatment and control group and the average number of bottles was quite a bit lower than the treatment group, even though it was still a bit above zero.

I decided, on a lark, to put the number of bottles received in as a covariate in my analysis. To my initial horror, the effect of the treatment disappeared when you adjusted for the number of bottles.

But I quickly realized that this, if anything, reinforced the conclusions of the study. The number of bottles given was part of the causal pathway. It occured after random assignment, and the intervention was deliberately designed to influence this intermediate variable. So the adjustment actually helped to explain why the intervention worked.

Adjusting for baseline measurements

Baseline = measurements prior to intervention
- Done to improve precision
- Can use baseline as a covariate
Change score is an alternative
- Also known as difference in differences (DID) model
- Possible regression to the mean

Many research studies include baseline measurements, measurement of the same outcome measure that you plan to use to compare the treatment and control groups. It often helps to make an assessment of the outcome BEFORE you implement any intervention. This is done mostly to improve precision, but it can also help control bias. If the intervention and control groups differ on the baseline measure, that is an indiciation that one group is more seriously ill at the start of the research than the other group.

You can consider the baseline value to be a covariate. Your outcome at baseline certainly does have some influence on your outcome at the end of the study. But there is no direct interest in the baseline measure itself.

A common alternative to using the baseline measure as a covariate is to compute change scores, the difference between the baseline measure of the outcome and the measure of the outcome at the end of the study. This is measuring the relative decline or improvement in health.

Use of change scores or difference in difference models is controversial. I like this approach but most of the research community is criticizes this choice. One criticism is that regression to the mean can possibly mess up the change score analysis.

Regression to the mean is the tendency for extremely low scores at one time point often are not quite as extreme when they are measured again, even if there is no change or intevention going on. Similarly extremely high scores at one time point often are not as extreme when measured again.

What did we learn in module 04?

Mathematical model, 1

Decompose \(\mu_{ij}\) into \(\mu + \alpha_i + \beta_j\)
- \(\alpha_i\) is the deviation for the ith level of first factor
- \(\beta_j\) is the deviation for the jth level of second factor
- Require \(\alpha_1=0\) and \(\beta_1=0\)
- \(\mu\) is the mean for the reference levels

Mathematical model, 2

\(Y_{ijk} = \mu + \alpha_i + \beta_j +\epsilon_{ijk}\)
- i=1,…,a levels of the first categorical variable
- j=1,…,b levels of the second categorical variable
- k=1,…,n replicates with first and second categories
Note: \(\mu, \alpha_i, \beta_j, \epsilon_{ijk}\) are population values

Mathematical model, 3

\(H_0:\ \alpha_i=0\) for all i

\(\ \)

\(H_0:\ \beta_j=0\) for all j

Parameter estimates for the two factor model

# A tibble: 14 × 5
   term        estimate std.error statistic p.value  
   <chr>          <dbl>     <dbl>     <dbl> <glue>   
 1 (Intercept)    4.72      1.50      3.14  p = 0.005
 2 MoonDuring     2.5       0.984     2.54  p = 0.019
 3 MoonAfter      0.542     0.984     0.550 p = 0.588
 4 MonthSep       4.03      1.97      2.05  p = 0.053
 5 MonthOct       3.73      1.97      1.90  p = 0.071
 6 MonthNov       3.70      1.97      1.88  p = 0.073
 7 MonthDec       2.63      1.97      1.34  p = 0.195
 8 MonthJan       5.03      1.97      2.56  p = 0.018
 9 MonthFeb       6.93      1.97      3.52  p = 0.002
10 MonthMar       8.57      1.97      4.35  p < 0.001
11 MonthApr      13.0       1.97      6.61  p < 0.001
12 MonthMay       8.60      1.97      4.37  p < 0.001
13 MonthJun       7.10      1.97      3.61  p = 0.002
14 MonthJul      11.0       1.97      5.61  p < 0.001

Analysis of variance table comparing the two factor model to the null model

# A tibble: 2 × 7
  term                     df.residual   rss    df sumsq statistic p.value  
  <chr>                          <dbl> <dbl> <dbl> <dbl>     <dbl> <glue>   
1 Admission ~ 1                     35  625.    NA   NA      NA    <NA>     
2 Admission ~ Moon + Month          22  128.    13  497.      6.58 p < 0.001

Analysis of variance table comparing the two factor model to the one factor model

# A tibble: 2 × 7
  term                     df.residual   rss    df sumsq statistic p.value  
  <chr>                          <dbl> <dbl> <dbl> <dbl>     <dbl> <glue>   
1 Admission ~ Moon                  33  583.    NA   NA      NA    <NA>     
2 Admission ~ Moon + Month          22  128.    11  456.      7.13 p < 0.001

R-squared values

# A tibble: 3 × 3
  model  r.squared deviance
  <glue>     <dbl>    <dbl>
1 m1        0          625.
2 m2        0.0664     583.
3 m3        0.795      128.

Tukey post hoc test

# A tibble: 3 × 7
  term  contrast      null.value estimate conf.low conf.high adj.p.value
  <chr> <chr>              <dbl>    <dbl>    <dbl>     <dbl> <chr>      
1 Moon  After-Before           0    0.542  -1.93        3.01 0.847      
2 Moon  During-Before          0    2.50    0.0280      4.97 0.047      
3 Moon  During-After           0    1.96   -0.514       4.43 0.138

What did we learn in module 05?

Mathematical model, 1

\(Y_{ijk}=\mu+\alpha_i+\beta_j+(\alpha \beta)_{ij}+\epsilon_{ijk}\)
- i=1,…,a j=1,…,b, k=1,…,n
If 1 is the reference category
- \(\alpha_1=0\)
- \(\beta_1=0\)
- \((\alpha \beta)_{1j}=0\)
- \((\alpha \beta)_{i1}=0\)

Mathematical model, 2

\(SS_A=\Sigma_i nb(\bar{Y}_{i..}-\bar{Y}_{...})^2\)
\(SS_B=\Sigma_i na(\bar{Y}_{.j.}-\bar{Y}_{...})^2\)
\(SS_{AB}=\Sigma_i \Sigma_j n(\bar{Y}_{ij.}-\bar{Y}_{i..}-\bar{Y}_{.j.}+ \bar{Y}_{...})^2\)
\(SS_E=\Sigma_i \Sigma_j \Sigma_k (Y_{ijk}-\bar{Y}_{ij.})^2\)
\(SS_T=\Sigma_i \Sigma_j \Sigma_k (Y_{ijk}-\bar{Y}_{...})^2\)

Test for an interaction

\(SS_{AB}\) has (a-1)(b-1) degrees of freedom
\(SS_E\) has ab(n-1) degrees of freedom
Accept \(H_0\) if \(F=\frac{MS_{AB}}{MS_E}\) is close to one
- In R, fit a model without an interaction
- Compare to a model with interaction
- Using the anova function

What did we learn in module 06?

Comparing two binary outcomes

Is there a difference in the proportion of deaths between male passengers and female passengers on the Titanic?
Is there difference in the proportion of patients finishing the full three doses of HPV vaccine between Black women and White women?
Does using a ng tube for feeding in pre-term infants increase the probability of successful breast feeding at six months?

Other comparisons involving a binary outcome

Is there are difference in the proportion of deaths between first class, second class, and third class passengers?
Does age influence the proportion of women finishing the full three doses of HPV vaccine?
Controlling for the mother’s age, does using a ng tube for feeding in pre-term infants increase the probability of successful breast feeding at six months?

Hypothesis framework

\(H_0:\ \pi_1=\pi_2\)
\(H_1:\ \pi_1=\pi_2\)
Compute \(\hat p_1\) and \(\hat p_2\) from samples
Accept \(H_0\) if \(\hat p_1-\hat p_2\) is close to zero.
- \(T=(\hat p_1-\hat p_2)/s.e.\)
- 95% CI: \((\hat p_1-\hat p_2) \pm Z_{\alpha/2}s.e.\)

The Titanic dataset

Rows: 1,313
Columns: 5
$ Name     <chr> "Allen, Miss Elisabeth Walton", "Allison, Miss Helen Loraine"…
$ PClass   <chr> "1st", "1st", "1st", "1st", "1st", "1st", "1st", "1st", "1st"…
$ Age      <dbl> 29.00, 2.00, 30.00, 25.00, 0.92, 47.00, 63.00, 39.00, 58.00, …
$ Sex      <chr> "female", "female", "male", "female", "male", "male", "female…
$ Survived <dbl> 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1…

Counts and percentages

        Survived
Sex      Yes  No
  female 308 154
  male   142 709

        Survived
Sex            Yes        No
  female 0.6666667 0.3333333
  male   0.1668625 0.8331375

Test for difference in proportions

# A tibble: 1 × 9
  estimate1 estimate2 statistic  p.value parameter conf.low conf.high method    
      <dbl>     <dbl>     <dbl>    <dbl>     <dbl>    <dbl>     <dbl> <chr>     
1     0.667     0.167      332. 3.43e-74         1    0.450     0.550 2-sample …
# ℹ 1 more variable: alternative <chr>

Chi-square test of independence, 1 of 2

Equivalent to test of two proportions
Lay out data in two by two table

\[\begin{matrix} & No\ event & Event \\ Treatment & O_{11} & O_{12}\\ Control & O_{21} & O_{22} \end{matrix}\]

Chi-square test of independence, 2 of 2

\[\begin{matrix} & No\ event & Event \\ Treatment & E_{11} = n_1 (1-\hat p_.) & E_{12}=n_1 \hat p_.\\ Control & E_{21} = n_2 (1-\hat p_.) & E_{22}=n_2 \hat p_. \end{matrix}\]

\(X^2=\Sigma \frac{(O_{ij}-E_{ij})^2}{E_{ij}}\)

Expected counts for Titanic

Observed counts

        Survived
Sex      Yes  No
  female 308 154
  male   142 709

Expected counts

        Survived
Sex           Yes       No
  female 158.3397 303.6603
  male   291.6603 559.3397

Chisquare test for Titanic

# A tibble: 1 × 4
  statistic  p.value parameter method                    
      <dbl>    <dbl>     <int> <chr>                     
1      332. 3.43e-74         1 Pearson's Chi-squared test

Odds ratio calculation

       No event  Event  Odds
Group1    a         b
Group2    c         d

Odds for group 1 = \(b/a\)
Odds for group 2 = \(d/c\)
Odds for group 1 = \(\frac{d/c}{b/a} = \frac{ad}{bc}\)
s.e.(log or) = \(\sqrt{\frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d}}\)

Titanic data

       Survived   Died  Total
Female   308      154     462
Male     142      709     851
Total    450      863   1,313

Titanic data, odds of death

       Survived   Died  Total  Odds
Female   308      154     462  2     to 1 against
Male     142      709     851  4.993 to 1 in favor
Total    450      863   1,313

Odds ratio = 4.993 / 0.5 = 9.986

Odds ratio for survival by sex

$data
        Survived
Sex      Yes  No Total
  female 308 154   462
  male   142 709   851
  Total  450 863  1313

$measure
        odds ratio with 95% C.I.
Sex      estimate    lower    upper
  female 1.000000       NA       NA
  male   9.956188 7.662525 13.00928

$p.value
        two-sided
Sex      midp.exact fisher.exact   chi.square
  female         NA           NA           NA
  male            0 4.826448e-74 3.425855e-74

$correction
[1] FALSE

attr(,"method")
[1] "median-unbiased estimate & mid-p exact CI"

What did we learn in module 07?

What is a diagnostic test?

Indication of disease
- Rapid, convenient, and/or inexpensive
Gold standard
- Indication of same disease
- Slow, inconvenient, and/or expensive

Examples of diagnostic tests, 1 of 2

Yale-Brown obsessive-compulsive scale
- Do you often feel sad or depressed?
SCOFF questionnaire
- Five yes/no questions
- Two or more yes responses

Some examples of diagnostic tests are:

The Yale-Brown obsessive-compulsive scale, a simple yes/no answer to the following question: Do you often feel sad or depressed? In a study of stroke patients at the Royal Liverpool and Broadgreen University Hospitals (BMJ 2001; 323: 1159), this test was shown to perform well compared to a more complex measure, the Montgomery Asberg depression rating scale.

The SCOFF questionnaire asks five yes/no questions to determine whether a patient has an eating disorder.

Do you ever make yourself sick because you feel uncomfortably full?
Do you worry you have lost control over how much you eat?
Have you recently lost more than one stone in a 3 month period?
Do you believe yourself to be fat when others say you are thin?
Would you say that food dominates your life?

Two or more yes answers is considered a positive test. In a study of 341 consecutive patients at two general practices in southwest London (BMJ 2002; 325: 755-756), these patients were given the SCOFF questionnaire and then a formal interview based on Diagnostic and Statistical Manual of Mental Disorders, (fourth edition). The interview lasted 10-15 minutes and the interviewer did not know that score on the SCOFF questionnaire. The SCOFF questionnaire produced results that were comparable to the formal interview.

More examples

Rectal bleeding as a sign of colorectal cancer
Electrocardiogram, QTc dispersion

What did we learn in module 08?

Survival analysis

Time to event models
- Death
- Relapse
- Rehospitalization
- Failure of medical device
- Pregnancy
Not every patient experiences the event
- These are censored observations

Survival analysis models are more properly called time to event models. You follow a group of patients from a certain time point and note the amount of time until they die.

Or the amount of time until they relapse. Or the amount of time until they need to be rehospitalized. Or the amount time until a medical device that you implanted in them fails.

I should note that while almost all of the events in a time to event model are bad, there are a few exceptions. In a study of couples with fertility problems, you might use a time to event model to study the time to pregnancy, a very happy event for any couple with fertility problems.

Mortality is the context under which time to event models were derived, so the term survival analysis has been used even when the event is different.

A key feature of survival analysis is that not every patient experiences the event. You should be glad that not everyone that you recruit for a clinical trial dies, but this adds a layer of complexity to the analysis.

The reasons for not experiencing the event can vary. Everyone dies, but maybe the event is death from cancer and if you patient gets hit by a bus, that patient does not experience the event. If the event is death from any cause, you still have to end the study in some time frame, and not everyone will die in that time frame. Not every patient gets rehospitalized, not every patient relapses, not every medical device fails.

A patient may drop out of a study, and you no longer are able to tell from that point onward whether that patient would have experienced the event sometime during the rest of your study.

When your patient does not die during the study, this is not a missing value. You have partial information. You know that the patient was alive for a certain amount of time. When you end the study within a certain time frame, you know that your patient dies at a time beyond the end date of your study. If your patient drops out after six months, you know that the patient survived for more than six months.

First fruit fly experiment, 1

data_dictionary: fly1.txt

description: |
  This dataset provides a simple example of what survival and censoring. It provides an inuitive explanation of estimation of survival probabilities.

vars:
  day:
    label: Time until death
    unit: days

First fruit fly experiment, 2

37, 40, 43, 44, 45, 47, 49, 54, 56, 58, 59, 60, 61, 62, 68, 70, 71, 72, 73, 75, 77, 79, 89, 94, 96

First fruit fly experiment, 3

First fruit fly experiment, 4

Second fruit fly experiment, 1

37, 40, 43, 44, 45, 47, 49, 54, 56, 58, 59, 60, 61, 62, 68, ??, ??, ??, ??, ??, ??, ??, ??, ??, ??

Second fruit fly experiment, 2

  day event
1  37     1
2  40     1
3  43     1
4  44     1
5  45     1
6  47     1
7  49     1
8  54     1
9  56     1

   day event
10  58     1
11  59     1
12  60     1
13  61     1
14  62     1
15  68     1
16  70     0
17  70     0
18  70     0

   day event
19  70     0
20  70     0
21  70     0
22  70     0
23  70     0
24  70     0
25  70     0

Second fruit fly experiment, 3

  day event   p
1  37     1 96%
2  40     1 92%
3  43     1 88%
4  44     1 84%
5  45     1 80%
6  47     1 76%
7  49     1 72%
8  54     1 68%
9  56     1 64%

   day event   p
10  58     1 60%
11  59     1 56%
12  60     1 52%
13  61     1 48%
14  62     1 44%
15  68     1 40%
16  70     0    
17  70     0    
18  70     0

   day event p
19  70     0  
20  70     0  
21  70     0  
22  70     0  
23  70     0  
24  70     0  
25  70     0

Second fruit fly experiment, 4

Third fruit fly experiment, 1

37, 40, 43, 44, 45, 47, 49, 54, 56, 58, 59, 60, 61, 62, 68, ??, 71, ??, ??, 75, ??, ??, 89, ??, 96

Third fruit fly experiment, 2

  day event
1  37     1
2  40     1
3  43     1
4  44     1
5  45     1
6  47     1
7  49     1
8  54     1
9  56     1

   day event
10  58     1
11  59     1
12  60     1
13  61     1
14  62     1
15  68     1
16  70     0
17  71     1
18  70     0

   day event
19  70     0
20  75     1
21  70     0
22  70     0
23  89     1
24  70     0
25  96     1

Third fruit fly experiment, 3

  day event   p
1  37     1 96%
2  40     1 92%
3  43     1 88%
4  44     1 84%
5  45     1 80%
6  47     1 76%
7  49     1 72%
8  54     1 68%
9  56     1 64%

   day event   p
10  58     1 60%
11  59     1 56%
12  60     1 52%
13  61     1 48%
14  62     1 44%
15  68     1 40%
16  70     0    
17  71     1 30%
18  70     0

   day event   p
19  70     0    
20  75     1 20%
21  70     0    
22  70     0    
23  89     1 10%
24  70     0    
25  96     1  0%

Third fruit fly experiment, 4

Here is a graph of the survival probability estimates from the third experiment.

These survival probabilities differ only slightly from the survival probabilities in the original experiment. This works out because the mechanism that caused us to lose information on six of the fruit flies was independent of their ultimate survival.

If the censoring mechanism were somehow related to survival prognosis, then you would have the possibility of serious bias in your estimates. Suppose for example, that only the toughest of flies (those with the most days left in their short lives) would have been able to escape. The flies destined to die on days 70, 71, 72, and 73, were already on their deathbeds and unable to fly at all, much les make a difficult escape. Then these censored values would not be randomly interspersed among the remaining survival times, but would constitute some of the larger values. But since these larger values would remain unobserved, you would underestimate survival probabilities beyond the 70th day.

This is known as informative censoring, and it happens more often that you might expect. Suppose someone drops out of a cancer mortality study because they are abandoning the drugs being studied in favor of laetrile treatments down in Mexico. Usually, this is a sign that the current drugs are not working well, so a censored observation here might represent a patient with a poorer prognosis. Excluding these patients would lead to an overestimate of survival probabilities.

Interpreting Kaplan-Meier plots, 1

Interpreting Kaplan-Meier plots, 2

Interpreting Kaplan-Meier plots, 3

Interpreting Kaplan-Meier plots, 4

What did we learn in module 09?

Meta-analysis

Quantitative pooling of results from multiple studies
- Multi-center study
  - Each center has a different protocol
  - Some centers do not share results
Contrast to systematic overview
- Careful review of multiple studies
- May or may not include quantitative pooling
Contrast to scoping review
- “Researchers may conduct scoping reviews instead of systematic reviews where the purpose of the review is to identify knowledge gaps, scope a body of literature, clarify concepts or to investigate research conduct.” Munn 2018

Meta-analysis is the quantitative pooling of results from multiple independently published research studies. I joke about how meta-analysis is a multi-center research study but with a couple of qualifications. First, each center gets to use a different protocol. Second, some centers do not share their results with you. This hints at a couple of important issues I will talk about in detail: heterogeneity and publication bias. When I describe it as a chaotic multi-center trial, it sounds terrible. Well, maybe, but we have still learned a lot from meta-analytic studies in a broad range of scientific and medical areas.

A more commonly used term is “systematic overview” which is a superset of meta-analysis. A systematic overview is systematic meaning that it uses a careful and transparently documented approach to identify all research studies associated with a particular issue. It may or may not include a quantitative analysis. Thus all meta-analyses are systematic reviews but not all systematic reviews are meta-analyses.

There’s a new term, scoping review which I am less familiar with. Here’s a nice quote from a paper by Zachary Munn et al in BMC Medical Research Methodology. It sounds almost like a scoping review is a systematic overview with the pre-specified intent to stop before any serious meta-analytic intents.

I am going to focus only on meta-analysis because I am a statistician and this is a statistics course, so we all love anything quantitative. But never forget that there is more to research than just its quantitative aspects.

Case study: Declining sperm counts

Meta-analysis published in 1992 in BMJ
- 62 studies
- 1938 through 1991
- 12 million/ml decline per decade
Criticism
- Early studies in North America, later studies more global
- Variations in collection methods, abstinence requirements
- Many competing analyses
Bad news/good news

In 1992, the British Medical Journal published a controversial meta-analysis. This study (BMJ 1992: 305(6854); 609-13) reviewed 61 papers published from 1938 and 1991 and showed that there was a significant decrease in sperm count and in seminal volume over this period of time. For example, a linear regression model on the pooled data provided an estimated average count of 113 million per ml in 1940 and 66 million per ml in 1990.

Several researchers (Fertil Steril 1996: 65(5); 1044-6 and Fertil Steril 1995: 63(4); 887-93) noted heterogeneity in this meta-analysis, a mixing of apples and oranges. Studies before 1970 were dominated by studies in the United States and particularly studies in New York. Studies after 1970 included many other locations including third world countries. Thus the early studies were United States apples. The later studies were international oranges. There was also substantial variation in collection methods, especially in the extent to which the subjects adhered to a minimum abstinence period.

The original meta-analysis and the criticisms of it highlight both the greatest weakness and the greatest strength of meta-analysis.

Meta-analysis is the quantitative pooling of data from studies with sometimes small and sometimes large disparities. Think of it as a multi-center trial where each center gets to use its own protocol and where some of the centers are left out.

On the other hand, a meta-analysis lays all the cards on the table. Sitting out in the open are all the methods for selecting studies, abstracting information, and combining the findings. Meta-analysis allows objective criticism of these overt methods and even allows replication of the research.

Contrast this to an invited editorial or commentary that provides a subjective summary of a research area. Even when the subjective summary is done well, you cannot effectively replicate the findings. Since a subjective review is a black box, the only way, it seems, to repudiate a subjective summary is to attack the messenger.

Major issues in meta-analysis

Heterogeneity
- Were apples combined with oranges?
Publication bias
- Were some apples left on the tree?
Study quality
- Were all the apples rotten?
Interpretability
- Did the pile of apples amount to more than just a hill of beans?

What did we learn in module 10?

Talk given to first year medical students

Topic also relevant to this class.
Only a few minor changes
- Different format for the “programming” assignment

Who am I?

Steve Simon

PhD Statistics, 1982, U Iowa
Teach in Biomedical and Health Informatics
- Previous jobs at CMH, CDC
Part-time independent statistical consultant (P.Mean Consulting)
Married to a Pediatric Cardiologist (retired)
Run 5K and 4 mile races

Obsessed with computers since 1972

Figure 1. Section on computer skills from my resume

Worked with health care applications since 1987

Recent positions
- Centers for Disease Control and Prevention (1987-1996)
- Children’s Mercy Hospital (1996-2008)
- UMKC School of Medicine (2008 to present)
But…
- I am not a doctor
- Still confused about many things
  - Example: Difference between good and bad cholesterol.

Quiz questions (1/3)

Why does Joel Best call statistics a social construct?

Statistics are misquoted often on social media.
Statistics are selected, shaped, and presented by human beings.
Statistics are used to promote socialism.
Statistics are dehumanizing.

Quiz questions (2/3)

What is the main philosophical foundation of empiricism?

Everything can be reduced to a mathematical equation.
Experiments can reveal the realities of the world.
Some questions are impossible to answer.
We construct our own reality based on our own lived experiences

Quiz questions (3/3)

What is a major problem with data science?

Data scientists rely on large amounts of data with uneven quality.
Models developed by data scientists can lead to loss of privacy.
Prediction models are a black box that can hide discriminatory intent.
All of the above.

First poll question

Figure 2. Quote from “Peggy Sue Got Married”

Second poll question

Figure 3. Images of various computers

Are Statisticians Gods?

I’m helping someone who wants an alternative statistical analysis to the one used by the principal investigator. I’m happy to help and will offer advice about why my approach may be better, but I was warned that the PI considers the analysis chosen to be ordained by the “Statistical Gods” at her place of work.

What did we learn in module 11?

Hierarchical data

Moving beyond the independence assumption
Correlation within clusters

Examples of hierarchical data, 1 of 2

Body parts
- Left eye/right eye
- Teeth
- Skin patches
Human families
Animal litters

A simple example of hierarchical data is when you select a group of patients and then make measurements on two or more parts of their body. You might, for example, put an eye drop medication in the left eye and a placebo drop in the right eye. You might apply different types of sealents on different teeth in a mouth. You might put different food allergens on different parts of a patient’s back.

You might select families from a population and make measurements on two or more members of the same family. Since family members share the same environment and have very similar genetics, any comparison made within a family is likely to be more precise.

Likewise, measurements on the animals from the same litter will be precise because of a shared inter-uterine environment prior to birth and shared feeding from the same mother before weaning.

Examples of hierarchical data, 2 of 2

Clinics/hospitals
Communities
Repeated measurements

Longitudinal data (topic for next module)

Measurements taken at different times
- Emphasis in changes over time

Between and within cluster comparisons

Positive correlation
- Improves precision of within cluster comparisons
- Hurts precision of between cluster comparisons
Example with litters
- Medication administered during pregnancy
- Medication administered after birth

Basic notation, 1 of 2

\(Y_{ij}\)
- i defines cluster
  - i=1,…,a
- j defines individual within cluster
  - j=1,…,n

Basic notation, 2 of 2

\(Y_{ij} = \mu+\alpha_i+\epsilon_{ij}\)
- \(\mu\) unknown constant
- \(\alpha_i\) is normally distributed
  - \(SD(\alpha_i)=\sigma_{between}\)
- \(\epsilon_{ij}\) is normally distributed
  - \(SD(\epsilon_i)=\sigma_{within}\)

Some basic results

\(SD(Y_{ij})=\sigma_{total} = \sqrt{\sigma^2_{between}+\sigma^2_{within}}\)
\(SD(\bar{Y}_{..})=\sqrt{\frac{\sigma^2_{between}}{a}+\frac{\sigma^2_{within}}{an}}\)
\(Corr(Y_{ij}, Y_{ik})=\frac{\sigma^2_{between}}{\sigma^2_{between}+\sigma^2_{within}}\)
- Intraclass correlation (ICC)

Expected mean squares, 1 of 2

\(MS(between) = \frac{1}{a-1}\Sigma n(\bar{Y}_{i.}-\bar{Y}_{..})^2\)
- \(E[MS(between)] = n\sigma^2_{between}+\sigma^2_{within}\)

Expected mean squares, 2 of 2

\(MS(within) = \frac{1}{a(n-1)}\Sigma\Sigma(\bar{Y}_{ij}-\bar{Y}_{i.})^2\)
- \(E[MS(within)] = \sigma^2_{within}\)

Variance components estimates

\(\hat\sigma_{between}^2=\frac{MS(between)-MS(within)}{n}\)
\(\hat\sigma_{within}^2=MS(within)\)

What did we learn in module 12?

Longitudinal data

Measurements taken at different times
- Emphasis in changes over time

In the previous module, I talked about hierarchical models and mentioned a particular case, longitudinal data, that I want to talk more about in this presentation.

Longitudinal data is similar to repeated measures data. With both, you measure the same subject repeatedly. With longitudinal data, often the emphasis is in changes that occur over time. Repeated measurements, in contrast, emphasize different treatments with the hope that the time gaps between the measurements are small enough that you don’t see changes over time.

The differences between longitudinal data, repeated measures data, or hierarchical data are subtle. Perhaps these are distinctions without a difference. I decided to separate out longitudinal data for a different module perhaps more out of the desire to split a complex topic into smaller bite-sized pieces.

Random intercepts model, 1

Simplest pattern for longitudinal data
\(Y_{ij},\ i=1,...,n;\ j=1,...,k\)
- n subjects, k time points
\(t_j\), time of jth measurement
- First time is often zero

Random intercepts model, 2

\(Y_{ij}=\beta_0+u_{0i}+\beta_1 t_j + \epsilon_{ij}\)
- \(\beta_0\) and \(\beta_1\) are unknown constants
- \(u_{0i}\) and \(\epsilon_{ij}\) are normally distributed
  - \(SD(u_{0i})=\sigma_{intercept}\)
  - \(SD(\epsilon_{ij})=\sigma_{error}\)

Random intercepts model, 3

\(SD(Y_{ij})=\sqrt{\sigma^2_{intercept}\ +\ \sigma^2_{error}}\)
\(Corr(Y_{ij}, Y_{im})=\frac{\sigma^2_{intercept}}{\sigma^2_{intercept}\ +\ \sigma^2_{error}}\)

Random intercepts illustrated, 1

Random intercepts illustrated, 2

Illustration of random intercepts with real data

What did we learn in module 13?

A simple example of Bayesian data analysis.

ECMO study
Treatment versus control, mortality endpoint
- Treatment: 28 of 29 babies survived
- Control: 6 of 10 babies survived
- Source: Jim Albert in the Journal of Statistics Education (1995, vol. 3 no. 3).

Bayesian data analysis seems hard, and it is. Even for me, I struggle with understanding Bayesian data analysis. In fairness, I must admit that much of my discomfort is just lack of experience with Bayesian methods. In fact, I’ve found that in some ways, Bayesian data analysis is simpler than classical data analysis. You, too, can understand Bayesian data analysis, even if you’ll never be an expert at it. There’s a wonderful example of Bayesian data analysis at work that is simple and fun. It’s taken directly from an article by Jim Albert in the Journal of Statistics Education (1995, vol. 3 no. 3) which is available on the web at www.amstat.org/publications/jse/v3n3/albert.html.

I want to use his second example, involving a comparison of ECMO to conventional therapy in the treatment of babies with severe respiratory failure. In this study, 28 of 29 babies assigned to ECMO survived and 6 of 10 babies assigned to conventional therapy survived. Refer to the Albert article for the source of the original data. Before I show how Jim Albert tackled a Bayesian analysis of this data, let me review the general paradigm of Bayesian data analysis.

Wikipedia introduction

P(H|E) = P(E|H) P(H) / P(E)
- H = hypothesis
- E = evidence
- P(H) = prior
- P(E|H) = likelihood
- P(H|E) = posterior

Prior distribution

Degree of belief
- Based on previous studies
- Subjective opinion (!?!)
Examples of subjective opinions
- Simpler is better
- Be cautious about subgroup analysis
- Biological mechanism adds evidence
Flat or non-informative prior

The first step is to specify P(H), which is called the prior probability. Specifying the prior probability is probably the one aspect of Bayesian data analysis that causes the most controversy. The prior probability represents the degree of belief that you have in a particular hypothesis prior to collection of your data. The prior distribution can incorporate data from previous related studies or it can incorporate subjective impressions of the researcher. What!?! you’re saying right now. Aren’t statistics supposed to remove the need for subjective opinions? There is a lot that can be written about this, but I would just like to note a few things.

First, it is impossible to totally remove subjective opinion from a data analysis. You can’t do research without adopting some informal rules. These rules may be reasonable, they may be supported to some extent by empirical data, but they are still applied in a largely subjective fashion. Here are some of the subjective beliefs that I use in my work:

You should always prefer a simple model to a complex model if both predict the data with the same level of precision.

You should be cautious about any subgroup finding that was not pre-specified in the research protocol.

if you can find a plausible biological mechanism, that adds credibility to your results.

Advocates of Bayesian data analysis will point out that use of prior distributions will force you to explicitly state some of the subjective opinions that you bring with you to the data analysis.

Second, the use of a range of prior distributions can help resolve controversies involving conflicting beliefs. For example, an important research question is whether a research finding should “close the book” to further research. If data indicates a negative result, and this result is negative even using an optimistic prior probability, then all researchers, even those with the most optimistic hopes for the therapy, should move on. Similarly, if the data indicates a positive result, and this result is positive even using a pessimistic prior probability, then it’s time for everyone to adopt the new therapy. Now, you shouldn’t let the research agenda be held hostage by extremely optimistic or pessimistic priors, but if any reasonable prior indicates the same final result, then any reasonable person should close the book on this research area.

Third, while Bayesian data analysis allows you to incorporate subjective opinions into your prior probability, it does not require you to incorporate subjectivity. Many Bayesian data analyses use what it called a diffuse or non-informative prior distribution. This is a prior distribution that is neither optimistic nor pessimistic, but spreads the probability more or less evenly across all hypotheses.